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Abstract. We consider the discontinuous Petrov-Galerkin (DPG) method, where the 

test space is normed by a modified graph norm. The modification scales one of the terms 

in the graph norm by an arbitrary positive scaling parameter. Studying the application 

O ! of the method to the Hclmholtz equation, we find that better results arc obtained, under 

•^T some circumstances, as the scaling parameter approaches a limiting value. We perform a 

dispersion analysis on the multiple interacting stencils that form the DPG method. The 
analysis shows that the discrete wavenumbers of the method are complex, explaining 
the numerically observed artificial dissipation in the computed wave approximations. 

^* Since the DPG method is a nonstandard least-squares Galerkin method, we compare its 






> 



% 



performance with a standard least-squares method. 



1. Introduction 



Discontinuous Petrov-Galerkin (DPG) methods were introduced in [8], [10]. The DPG 
methods minimize a residual norm, so they belong to the class of least-squares Galerkin 
methods [31 [JJ [H], although the functional setting in DPG methods is nonstandard. In 
this paper, we introduce an arbitrary parameter e > into the definition of the norm in 
which the residual is minimized. We study the properties of the resulting family of DPG 



ON 

r methods when applied to the Helmholtz equation. 

The DPG framework has already been applied to the Helmholtz equation in [12] . An 
error analysis with optimal error estimates was presented there. There are two major 
differences in the content of this paper and [12] • The first is the introduction of the above 
mentioned parameter, e. When e = 1, the method here reduces to that in [T2J. The use of 

^ such scaling parameters was already advocated in [TT] based on numerical experience. In 

this paper, we shall provide a theoretical basis for its use. The second major difference 
with [12j is that in this contribution we perform a dispersion analysis of the DPG method 
with the e scaling. We thus discover several important properties of the method as e is 
varied. 

Least-squares Galerkin methods are popular methods in scientific computation [3j [17] . 
They yield Hermitian positive definite systems, notwithstanding the indefiniteness of the 
underlying problem. Hence they are attractive from the point of view of solver design 
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and many works have focused on this subject [T5j IT5] . However, as we shall shortly 
see in detail, for wave propagation problems, they yield solutions with heavy artificial 
dissipation. Since the DPG method is of the least-squares type, it also suffers from this 
problem. One of the goals of this paper is to show that by means of the ^-scaling, we can 
rectify this problem to some extent. 

To explain this issue further, let us fix the specific boundary value problem we shall 
consider. Let A : #(div, Q) x H X (Q) -> I?{Q) N x L 2 {Q) denote the Helmholtz wave 
operator defined by 

(1) A(v,r)) - (uov + Vr),iu)r) + V -v). 

Here i denotes the imaginary unit, u is the wavenumber, and Q is a bounded open 
connected domain with Lipschitz boundary. All function spaces in this paper are over 
the complex field C The Helmholtz equation takes the form A(u,<j)) = f, for some f € 
L 2 (f2) N x L 2 (Q). Although, we consider a general f in this paper, in typical applications, 
f = (0, /) with / € L 2 (f2), in which case, eliminating the vector component u , we recover 
the usual second order form of the Helmholtz equation, 

-A(f> - u 2 cf) = iuf, on Q. 

This must be supplemented with boundary conditions. The DPG method for the case 
of the impedance boundary conditions iojcj) + dtfi/dn = on dQ was discussed in [12] , but 
other boundary conditions are equally well admissible. In the present work, we consider 
the Dirichlet boundary condition 

(2) = 0, on dfl. 

To deal with this boundary condition, we will need the space 

(3) R = H(div,Q)xH^(Q), 
Thus, our boundary value problem reads as follows: 

(4) Find (u,(j)) e R satisfying A(u , (/)) = f. 

It is well known |TB] that except for u in E, an isolated countable set of real values, this 
problem has a unique solution. We assume henceforth that u is not in E. 

Before studying the DPG method for Q, it is instructive to examine the simpler 1? 
least-squares Galerkin method. Set Rh c R to the Cartesian product of the lowest order 
Raviart-Thomas and Lagrange spaces, together with the boundary condition in R. The 
method finds (u^,0^) e Rh such that 

(5) (u i s , (f)l) = arg mm || f - Aw || . 



wei? 



h 



Throughout, | • | denotes the L 2 (i7) norm, or the natural norm in the Cartesian product 
of several L 2 (i7) component spaces. The method (|5| belongs to the so-called FOSLS [7] 
class of methods. 

Although ^ appears at first sight to be a reasonable method, computations yield solu- 
tions with artificial dissipation. For example, suppose we use (|5]), appropriately modified 
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to include nonhomogeneous boundary conditions, to approximate a plane wave propagat- 
ing at angle 9 = n/8 in the unit square. A comparison between the real parts of the exact 



solution (in Figure la) and the computed solution (in Figure lb) shows that the computed 



solution dissipates at interior mesh points. The same behavior is observed for the lowest 



order DPG method with e - 1 in Figure lc (see £2.4 for the definition of r therein and 



Section H for a full discussion of the lowest order DPG method). The same method with 



e = 10~ 6 however gave a solution (in Figure Id) that is visually indistinguishable from 
the exact solution. Note that, for the DPG method with e = 1, the numerical results 
presented in [TJ] show much better performance, because slightly higher order spaces 
were used there. Instead, in this paper, we have chosen to study the DPG method with 
the lowest possible order of approximation spaces to reveal the essential difficulties with 
minimal computational effort. 



The situation in Figures lb and lc improves when more elements per wavelength are 
used. This is not surprising in view of the asymptotic error estimates of the methods. To 
give an example of such an error estimate, consider the case of the impedance boundary 
conditions considered in [12J. It is proven there that there is a constant C > 0, independent 
of u and mesh size h, such that the lowest order DPG solution (uh,<f>h) satisfies 

(6) \\u -u h \\ + \\(f)-<f) h \\ <Cu 2 h 

for a plane wave solution. A critical ingredient in this analysis is the estimate 

(7) ||m/|1<C"|Aw|, 

which, as shown in [12j Lemmas 4.2 and 4.3], holds for all w in the analogue of R with 
impedance boundary conditions. Although the analysis in [12] was for the impedance 
boundary condition, similar techniques apply to the Dirichlet boundary condition as well, 
leading to (J6J). As more elements per wavelength are used, uih decreases, so (|6| guarantees 



that the situation in Figure [Tcj will improve. 

The analysis for the L 2 least-squares method is easier than the above-mentioned DPG 
analysis. Indeed, by ([5}, \\f - A(v,h,<f%)\\ < \\A(u - w h ,<p - ip h )\\ for any (w h ,ip h ) e R h . 
Hence, applying (J7|) to the error e = (u - u\,cj) - 0^) and noting that the residual is 
Ae = f- A(u^,(fi l £), we obtain |e| < C"||A(w -Wh,4>-iph)\- By standard approximation 
estimates, we then conclude that there is a C > independent of ui and h such that 

(8) \\u~ut\\ + H-$\\<Cu 2 h. 

This simple technique of analysis of L 2 -based least-squares methods is well-known (see e.g., 
[T?| pp. 70-71]). As with (Jsl) , the estimate pi) implies that as the number of elements per 



wavelength is increased, uih decreases, and the situation in Figure lb must improve. 



Yet, Figures lb and lc show that these methods fail to be competitive with standard 
methods in accuracy for small number of elements per wavelength. The figures also 
illustrate one of the difficulties with asymptotic error estimates like (J7l) and (J8J). Having 
little knowledge of the size of C, we cannot predict the performance of the method on 



coarse meshes. Motivated by this difficulty, one of the theorems we present (Theorem 3.1 ) 



will give a better idea of the constant involved as e -*■ 0. Also note that the above 
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(a) A plane wave propagating at angle ir/8 





(b) L 2 least-squares solution 




(c) Numerical traces from the lowest order (d) Numerical traces from the lowest order 
DPG method with e - 1 and r - 3 DPG method with e - 1CT 6 and r - 3 

Figure 1. Approximations to a plane wave computed using a uniform 
mesh of square elements of size h = 1/48 (about sixteen elements per wave- 



length). Artificial dissipation is visible in Figures lb and lc 



indicated error analyses does not give us a quantitative measure of differences in wave 
speeds between the computed and exact waves. This motivates the dispersion analysis we 
present in this paper, which will address the issue of wave speed discrepancies. 

We should note that there are alternative methods of the least-squares type that exhibit 
better performance than the standard L 2 -based least squares method. Some are based 
on adding further terms to the residual to be minimized (e.g., to control the curl of 
the vector equation [18]). Another avenue explored by others, and closer to the subject 
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of this paper, is the idea of minimizing the residual in a dual norm jH |S|. The main 
difference with our method is that our dual norms are locally computable in contrast to 
their nonlocal norms. This is achieved by using an ultraweak variational setting. The 
domain and codomain of the operator in the least-squares minimization associated to the 
DPG method are nonstandard, as we shall see next. 

2. The DPG method for the Helmholtz equation 

In this section, we briefly review the method for the Helmholtz equation introduced 
in [T2]. We then show exactly where the parameter e is introduced to get the variant of 
the method that we intend to study. 

Let flh be a disjoint partitioning of Q into open elements K such that Q = u K( ,Q h K. 
The shape of the mesh elements in i?^ is unimportant for now, except that we require 
their boundaries dK to be Lipschitz so that traces make sense. Let 

(9) V = H(div,f2 h )xH 1 ((2 h ), 
where 

H(div,Q h ) = {f: f\ K eH(div,K), VK*Q h }, 

H\Q h ) = {v : v\ K € H\K), VK € Q h ). 

Let Ah ■ V -> L 2 (Q) N x L 2 (i?) be defined in the same way as A in (pj), except the 
derivatives are taken element by element, i.e., on each K e i?^, we have Ah(v,r])\K = 
(iuv\ K + Vtj\k, 10JT]\k + V • v\ K ). 

2.1. Integration by parts. The following basic formula that we shall use is obtained 
simply by integrating by parts each of the derivatives involved: 

(10) / A(w ,ip) ■ (v,r)) - - / (w , ■0) • A(v, rj) + / (w-n)r}+ ip(v-n), 

JD JD JdD ' JdD 

for smooth functions (w ,ip) and (v,rj) and domains D with Lipschitz boundary. Above, 
overlines denote complex conjugations and the integrals use the appropriate Lebesgue 
measure. Note that we use the notation n throughout to generically denote the outward 
unit normal on various domains - the specific domain will be clear from context - e.g., 
in (10), it is D. Introducing the following abbreviated notations for tuples w = (w ,ip) 
and v = (v,T]), 

(w,v) h = Y, / w -v + *Pv, 

KeH h JK 

((w,v)) h = V / (w-n)f)+ ip(v-n), 

(M, JdK JdK 



KeQ h 



we can rewrite (10), applied element by element, as 



(11) (Aw,v) h = -(w,A h v) h + ((w,v)) h . 



By density, (11) holds for all w e H(div 1 Q) x H l (Q) and all v e V. Then, ((•,•)) ^ must 



be interpreted using the appropriate duality pairing as the last term in (11) contains 
interelement traces on dQh = {dK : K e Qh)- 
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It will be convenient to introduce notation for such traces: Define 

tr h : if(div, Q) x H\n) -> n H- x l 2 (dK)n x H l ' 2 {dK) 

K 

as follows. For any (w ,ip) e H(div,f2) x if 1 (fi), the restriction of tr h(w ,ip) on the 
boundary of any mesh element OK takes the form ((w ■ n)n\dK,ip\dK) e H^ 1 / 2 (dK)n x 
H l ' 2 (dK). Although the meaning of H~ 1 / 2 (dK)n is more or less self-evident, to include 
a proper definition, first let Z denote the space of all functions of the form £n where £ is 
in if 1 / 2 (9i^), normed by ||£n|z = 1 £ || h V 2 (8k) ■ Let Z' denote the dual space of Z. Now, 
consider the map Mq - (q -n)n\dK, defined for smooth functions q on K. Since 



f Mq-in= f (q-n)£ 

JdK JdK 



(the left and right hand sides extend to duality pairings in Z and H x / 2 (dK), respectively), 
the standard trace theory implies that M can be extended to a continuous linear operator 
M : H (div, K) -> Z'. The range of M is what we denoted by ' L R- x l 2 {dK)nP Throughout 
this paper, functions in H^ 1 / 2 (dK)n appear together with a dot product with n, so we 
could equally well consider the standard space H~^ 2 (dK), but the notation simplifies 
with the former. In particular, with this notation, tr^w ,ip) is a single-valued function 
on the element interfaces since (w ,ip) is globally in H(div, Q) x iJ 1 (fi). 

2.2. An ultraweak formulation. The boundary value problem we wish to approximate 
is (13). Recall the definition of R in d3|. To deal with the Dirichlet boundary condition, 



we will need the trace space 

(12) Q = ti h (R). 

To derive the DPG method for 

(13a) A(u ,</>) = f, on i?, 

(13b) = 0, on dil, 

we use the integration parts by formula ( [TTj ) to get 

-((« , (f)),A h (v, rj)) h + {{tr h (u , 0), (v, 7]))) h = (f, (v, 7])} h 

for all (v,f]) € V. Now we let the trace tih(u,(j)) be an independent unknown (u,4>) in Q. 
Defining the bilinear form b((u,4>,u,4>),(v,r})) = -((u,<f)),A h (v,r)))h+ ((u,4>),(v,r)))) h , 
we obtain the ultraweak formulation of [12] : Find u = (u ,4>,u, <p) in 

U = L 2 {Q) N xL 2 (Q)xQ 

satisfying 

(14) b(u,v) = (f,v) h , VvzV. 

The wellposedness of this formulation was proved in [12] for the case of impedance bound- 
ary conditions. As is customary, we refer to the solution component u as the numerical 
flux and as the numerical trace. 
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2.3. The e-DPG method. Let U^ c U be a finite dimensional trial space. The DPG 
method finds Uh in Uh satisfying 

(15) b(u h ,v h ) = (f,v h ) h , 
for all v h in the test space Vh, defined by 

(16) V h = TU h , 
where T : U -*■ V is defined by 

(17) (Tw,v) v = b(w,v), VveV, 

and the V-inner product (■, -)y is the inner product generated by the norm 

n s"i L/ii 2 - ii /i ■/ii 2 j. C-2IU/II 2 



Here e > is an arbitrary scaling parameter. Note that when s = 1, (18) defines a graph 
norm on V. The case e = 1, analyzed in [12] . is the standard DPG method. In the next 
section, we will adapt the analysis of [32] to the case of the variable e, which we refer to 
as the "e-DPG method." 

It is easy to reformulate the e-DPG method as a residual minimization problem. (All 



DPG methods with test spaces as in (17) minimize a residual as already pointed out 
in [IH]-) Letting V denote the dual space of V, normed with | • \\v, we define F e V 
by F(v) = (f, v) h . Then letting B :U -*-V denote the operator generated by the above- 
defined &(•, •), i.e., Bw(v) = b(w, v) for all w e U and v e V, one can immediately see that 



Uh solves ( 15 ) if and only if 



Uh = arg min \\Bwh - F\\y. 

This norm highlights the difference between the DPG method and the previously discussed 
standard L 2 -based least-squares method ^. 



2.4. Inexactly computed test spaces. A basis for the test space Vh, defined in (16 
can be obtained by applying T to a basis of Uh- One application of T requires solving (17) 
which although local (calculable element by element), is still an infinite dimensional prob- 
lem. Accordingly a practical version of the e-DPG method uses a finite dimensional 
subspace V r c V and replaces T by T r : U -*■ V r defined by 

(19) (T r w,v) v = b(w,v), VveV r . 

In computations, we then use, in place of Vh, the inexactly computed test space V£ = T r Uh, 
i.e., the practical DPG method finds u r h in Uh satisfying 

(20) b(u h ,v) = {f,v) h , We^. 

For the Helmholtz example, we set V r as follows: Let Q/ m denote the space of polynomials 
of degree at most I and m in X\ and x 2 , resp. Let RT r = Q r>r -\ x Q r -i, r denote the Raviart- 
Thomas subspace of H(div, K). We set 

V r = {v : v\ K e RT r x Q rjr }. 
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Clearly, V r £ H (div, Q h ) x H 1 (i7 h ). Later, we shall solve (20) using r > 2 and report the 



numerical results. It is easy to see using the Fortin operators developed in [T5] that T r 



is injective for r > 2, which implies that (20) yields a positive definite system. However 



a complete analysis using [15] tracking uj and r dependencies, remains to be developed, 
and is not the subject of this paper. 

3. Analysis of the e-DPG method 
The purpose of this section is to study how the stability constant of the e-DPG 



method (15) depends on e. The analysis in this section provides the theoretical moti- 



vation to introduce the scaling by e into the DPG setting. 

3.1. Assumption. The analysis is under the already placed assumption that the bound- 



ary value problem (13) is uniquely solvable. We will now need a quantitative form of this 



assumption. Namely, there is a constant C(uo) > 0, possibly depending on u, such that 



the solution of (13) satisfies 

\\(u,<j>)\\<C(u)\\f\\. 

One expects C(co) to become large as u approaches any of the resonances in S. For any 
(r ,if)) 6 R, choosing f = A(r ,ip) and applying the above inequality, we obtain 

(21) \\(f^)\\<C(u)\\A(r^)\\, V(f,^)€ J R. 

This is the form in which we will use the assumption. 

Note that in the case of the impedance boundary condition, the unique solvability as- 
sumption can be easily verified [20J for all u. Furthermore, when that boundary condition 



is imposed, for instance, on the boundary of a convex domain, the estimate (21) is proved 
in [12], Lemmas 4.2 and 4.3] using a result of [20J . The resulting constant C{uj) is bounded 
independently of uj. However, we cannot expect this independence to hold for the Dirichlet 
boundary condition ^ we are presently considering. 

Finally, let us note that the ensuing analysis applies equally well to the impedance 
boundary condition: We only need to replace the space R considered here by that in [12] 



and assume (21) for all functions in the revised R. 



3.2. Quasioptimality. It is well-known that if there are positive constants C\ and Ci 
such that 

(22) CVHk<sup „ I, <C 2 \\v\\v, VveV, 

well \\ W \\U 



then a quasioptimal error estimate 

C 2 
(23) \\u- Uhlu < — inf |i7-w/||t/ 

d weU h 

holds. This follows from [12j Theorem 2.1], or from the more general result of [151 Theo- 
rem 2.1], after noting that the following uniqueness condition holds: Any w e U satisfying 
b(w,v) = for all v e V vanishes. (Since this uniqueness condition can be proved as 
in [T2l Lemma 4.1], we shall not dwell on it here.) 
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Accordingly, the remainder of this section is devoted to proving (22), tracking the 



dependence of constants with e, and using the ?7-norm we define below. First, let 

\\(f,ij)\\ R =-\\A(r^)\\. 



By virtue of (21), this is clearly a norm under which the space R, defined in (pi), is 
complete. The space Q in (12) is normed by the quotient norm, i.e., for any q eQ, 



\\cj\\q =inf {|/"|fl : for all r e R such that tr^r = q} . 

The function in R which achieves the infimum above defines an extension operator E : 
Q -> R that is a continuous right inverse of tr^ and satisfies 

/ c\ a\ II i — i / ^ II II ^ II 

(24) \\Eq\ R =\\q\\ Q . 

With these notations, we can now define the norm on the trial space by 

\{w,^,w,^)f u = \(w,^)f + \\(w, 



The following theorem is proved by extending the ideas in [12] to the £-DPG method 
Theorem 3.1. Suppose @ holds and let c = C(u) (c(w)e/2 + y/l + C(u/) 2 e 2 / 4 )- Th 



en 



the inf-sup condition in (22) holds with G\ = 1/vTTce and the continuity condition 
in (22) holds with C% = yl + ce. Hence, the DPG solution admits the error estimate 



\u- Uh\u <{! + ce) inf |u-w|t/. 

wsU h 

Proof. We first prove the continuity estimate. Let (w, q) e U and let v € V. We use the 



abbreviated notations q = (w,ip), w = (w,ip), and v = (v,rj). By (21) and (24), 

(25) \\Eq\\<C(u)£\\q\\ Q , \\AEq\\ = e\\q\\ Q . 

The extension E can be used to rewrite b((w, q),v) = -(w, A h v) h +(Eq ,A h v) h +(AEq , v) h . 



Then, applying the Cauchy-Schwarz inequality, and using (25), we have 
\b((w,q),v)\ < \\w\\\\A h v\\ +C(uj)e\\q\\ Q \\A h v\\ + e\\q \\ Q \\v\\ 

(26) <(M 2 + \\q\\l) 1/2 t, 

where t 2 = \\A h v\\ 2 + (C(u)e\\A h v\\ + e\\v\\) . With a = C(oj)e\\A h v\\ and b = e\\v\\ we apply 
the inequality (a + b) 2 < (l + a 2 )a 2 + (l + a~ 2 )b 2 to obtain 

t 2 < (1 + (1 + a 2 )C(u) 2 e 2 ) \\A h v\\ 2 + (1 + «- 2 ) £ 2 |>/|| 2 , 

for any a > 0. Setting a 2 = -1/2 + \/l/4 + C(a;)" 2 e" 2 , so that 

(27) (l + a 2 )C(a;)V = a- 2 = c£ 



with c as in the statement of the theorem. Hence, t 2 < (1 + ce)|v|y. Returning to (26) 

\b((w,q),v)\ <C 2 \\(w J q)\\ u \\v\\ v . 



with Ci = vl + ce. This verifies the upper inequality of (22). 
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To prove the lower inequality of (22), let r be the unique function in R satisfying Ar - v 



for any given v € V. Then, by (21 ), 

(28) 



<C(u) 



Also, since \\Ar\\ = \\v\\, letting f = tr^r, we have 



(29) 



\Q 



-\\AEf\\ < -\\Ar\ 
e e 



By (11), we have (Ar,v) h = -(r,A h v) h + {(f,v)) h , so 



£ 2 ||i/|| 2 + \\A h v\ 



e 2 b((z,r),v), 



(30) \\v | 

where z = r - e^A^v, a function that can be bounded using p8J ), as follows: 

||z|| 2 < (1 + a 2 )\\r | 2 + (1 + cT 2 )£- 4 || A h v\\ 2 

< (1 + a 2 )C(a;) 2 ||i/|| 2 + (1 + oT 2 K 4 || A>/|| 2 , 



for any a > 0. Choosing a as in ( |27[ ) and using (28)-(29), 

12 



l(z,r)| 



c/ 



+ £ \\ r \\Q 



(31) 



<(l + (l + « 2 )C(u;)V)£ 2 ||i/|| 2 + (l + a- 2 )||A h i/| 

<(l + C£)(£ 2 |l/| 2 +|^l/| 2 ). 



Returning to ( 30 ) , we now have 
II. .1,2 K( z ^),v) 



( z ) r )|c/ \«I/ |X|| 17 / 



by virtue of (31), verifying the lower inequality of (22) with C\ = ljyl + ce. 



D 



Remark 3.2. Although we presented the above result only for the Helmholtz equation, 
the ideas apply more generally. It seems possible to prove a similar result abstractly, e.g., 
using the abstract setting in [6], for any DPG application that uses a scaled graph norm 
analogous to (18) (with the wave operator Ah replaced by suitable others). 



3.3. Discussion. Theorem 3A shows that the use of the e-scaling in the test norm can 
ameliorate some stability problems, e.g., those that can arise from large C(cu). 

Observe that the best possible value for the constant Ci\C\ in (23) is 1. Indeed, if 
C2/C1 equals 1, then the computed solution u^ coincides with the best approximation to 
u from Uh- Theorem 3.1 shows that the quasioptimality constant of the DPG method 
approaches the ideal value of 1 as e -> 0. 

However, since the norms depend on e, we must further examine the components of the 
error separately, by defining 



(32a) 
(32b) 



\u -u h \\ 2 + 



e 2 = \\AE(u-u h ,<J)- ( J) h )\\ 2 . 
The estimate of Theorem |3.1| implies that 



(33) 



e + 



<(l + ce) 2 (a 2 + ^) 
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where a and a are the best approximation errors defined by 

(34) a 2 = inf \\u -w\\ 2 + 1|0- V|| 2 , 

(w,ip,0,0)eU h 

a 2 = inf |AB(u-t&,0-^)| 2 . 

(0,0,u),^)eC/ h 

Note that E is independent of e. 

We want to compare the error bounds for the numerical fluxes and traces in the e = 1 
case with the case of < e « 1. To distinguish these cases we will denote the error defined 



in (32b) by e\ when e = 1. Clearly, (33) implies 



(35) ej < (1 + cf (a 2 + a 2 ) . 



For the other case, (33) implies, after multiplying through by e 2 , 

e 2 < (l + ce) 2 (e 2 a 2 + a 2 ). 



Comparing this with (35), and noting that a and a remain the same for different e, we find 
that the DPG errors for fluxes and traces admit a better bound for smaller e. Whether the 
actually observed numerical error improves, will be investigated through the dispersion 
analysis presented in a later section, as well as in the next subsection. 



3.4. Numerical illustration. Theorem 3.1 partially explains a numerical observation 



we now report. We implemented the e-DPG method by setting the parameter r = 3 (see 



2.4) and computed u r h - (u r h , cj) r h , u r h: 4> r h ) . In analogy with (32), define the discretization 



errors e r and e r by e 2 = \\u - u r h \ 2 + ||0 - <jf h \ 2 and e 2 = \\AE(u - u r hl (f) - 4> r h ) || 2 . Although 

1/2 



Theorem |3.1| suggests an investigation of 

I"" u hlu (e 2 + (e r /ey 



inf weC/h fiy-i/i/lc/ \a 2 + (a/e) 2 ) 

due to the difficulty of applying the extension operator E in practice, we have investigated 

the ratio e r /a as a function of u. Recall that a is the L 2 (i7) best approximation error 

defined in p4| ), so e r /a measures how close the discretization errors are to the best possible. 

For a range of wavenumbers u, we chose the data f = (0, /) so that the exact solution 



to (13) on the unit square would be (u,(j)) = (-^V</>, 0), with <p = x(l - x)y(l - y). 
Each resulting boundary value problem was then solved using the e-DPG method with 
e = 10~ n ,n = 0, 1,2, 3, 4, on a fixed mesh of h = 1/16 and the corresponding discretization 
errors e r were collected. 

The resulting ratios e r /a are plotted as a function of u in Figure [2] for a few e values. 
First of all, observe that the graph of the ratio begins close to the optimal value of one 
for all e values in the figure. Next, observe that the ratio spikes up as u approaches the 
exact resonance value u> = 7r%/2 « 4.44, where C(ui) is infinity. It is interesting to look 
at the points near (but not at) the resonance. Observe that as e is decreased, the DPG 
method exhibits a "regularizing" effect at points near the resonance: E.g., at oj = 5, the 
values of e r /a are closer to 1 for smaller e. It therefore seems advantageous to use smaller 
e for problems near resonance. 
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Figure 2. The regularizing effect of e-DPG method as seen from a plot of 
the ratio e r /a near a resonance. 



The theoretical explanation for this numerical observation would be complete (by virtue 
of Theorem 3.1 ), if we had computed using the exact DPG test spaces (r = oo), instead of 
the inexactly computed spaces (r = 3). Certain discrete effects arising due to this inexact 
computation of test spaces will be presented in a later section. 



4. Lowest order stencil 

We now consider the example of square two-dimensional elements. The lowest order 
case of the DPG method is obtained using Q(dK) = {(w, , 0) : w is constant on each edge of 
dK, -if) is linear on each edge of dK, and ip is continuous on dK }. Let S(K) = {(w, ip) ■ w 
and ijj are constants (vector and scalar, resp.) functions on K}. We consider the DPG 
method (with e) using the lowest order global trial space 

U h = S h x Q hj 

where Qh = {f € Q ■ r\g K e Q(dK) for all mesh elements K} and Sh = {w ■ w\ K e S(K) for 
all mesh elements K}. 

Let x e denote the indicator function of an edge e. If a denotes a vertex of the square 
element K, let <p a denote the bilinear function that equals one at a and equals zero 
at the other three vertices of K. Let <p a = <p a \dK- The collection of eight functions of 
the form (0,(f> a ) and (x e ,0), one for each vertex, and one for each edge of K, forms a 
basis for Q(dK). We distinguish between the horizontal and vertical edges, because the 
unknowns there approximate different components of the velocity u . Accordingly, we will 
denote by Xe the indicator function of a horizontal edge and by xl the indicator function 
of a vertical edge. 
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(a) 21-point stencil 



(b) 13-point stencil 
Figure 3. Stencils 



(c) 13-point stencil 



The local 11 x 11 element DPG matrix is denned using a basis for Q(dK) and S(K). 
(While a basis for Q(dK) is obtained as mentioned above, a basis for S(K) is trivially 
obtained by three indicator functions.) If we enumerate the 11 basis functions as e;, 
% - 1, ... 11, then the element DPG matrix is defined by 



(36) 



Bij = b{e j ,T r e i ) 



where T r is as defined in (19). Since this matrix depends on u and e, we will write 



B = B(oj,e). In our computations, we do not use any specialized basis for V r to compute 
the action of T r , so to overcome round-off problems due to ill-conditioned local matrices, 
we resorted to high precision arithmetic for these local computations. 

To show how B can be computed by mapping, let K = [0,1] 2 . For any square K 
of side length h, there is a translation vector bx such that the K - bx = hK. For any 
(scalar or vector) function v on K, let v on K be the mapped function obtained by 



v(x) = v(hx + bx)- Let us denote the matrix computed using (36), but using the mapped 
basis functions e\ on K, by B(u,e). Then by a change of variables, it is easy to see that 



(37) 



B(u,e) = h 2 B(coh,eh). 



Thus we may compute local DPG matrix by scaling the DPG matrix on the fixed reference 
element K obtained using the normalized wavenumber uh and scaling parameter eh. It is 
enough to compute the element matrix B using high precision arithmetic for the ensuing 
dispersion analysis. 

Next, we eliminate the three interior variables of S(K) and consider the condensed 
8x8 element stiffness matrix for the variables in Q(dK). At this stage it will be useful to 
classify these eight variables (unknowns) into three categories: (1) Unknowns at vertices a 
(which are the coefficients multiplying the basis function a ) denoted by "•" , (2) unknowns 
on horizontal edges (coefficients multiplying Xe) denoted by "1"', and (3) unknowns on 
vertical edges (coefficients multiplying the corresponding x v e ) denoted by "_»" . The normal 
vectors on all horizontal and vertical edges are fixed to be (0, 1) and (1,0), respectively, 
corresponding to the direction of the above-indicated arrows. 

Now suppose the mesh is a uniform mesh of congruent square elements. Assembling 
the above-described condensed 8x8 element matrices on such a mesh, we obtain a global 
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system where the interior variables are all condensed out. The resulting equations can be 
represented using the stencils in Figure [3j A row of the matrix system corresponding to 
an unknown of the type "•" , connects to unknowns of the same type at other vertices, as 



well as unknowns of the other two types, as shown in the 21-point stencil in Figure 3a 



Similarly, the unknowns of the type "1*" and "_>" connect to other unknowns in the 13- 



point stencils depicted in Figures 3b and 3c, respectively. These stencils will form the 



basis of our dispersion analysis next. 



5. Dispersion analysis 

This section is devoted to a numerical study of the DPG method with e, by means of 
a dispersion analysis. The dispersion analysis is motivated by [13]. Details on dispersion 
analyses applied to standard methods can be found in [1] and the extensive bibliography 
presented therein. 

5.1. The approach. To briefly adapt the approach of [13] to fit our context, we consider 
a general method for the homogeneous Helmholtz equation on an infinite uniform lattice 
(/iZ) 2 . Suppose the method has S different types of nodes on this lattice, some falling 
in between the lattice points, each corresponding to a different type of variable, with its 
own stencil (and hence its own equation). All nodes of the t th type (t = l,2,...,S) are 
assumed to be of the form jh where j lies in an infinite subset of (Z/2) 2 . The solution 
value at a general node jh of the t th type is denoted by ip t j- Note that methods with 
multiple solution components are accommodated using the above mentioned node types. 

The t th stencil, centered around jh, consists of a finite number of nodes, some of which 
belong to the t th stencil, and the remaining belong to other stencils. Suppose we have finite 
index sets J s c (Z/2) 2 , for each s = 1,2, ... ,S, such that all the nodes of the t th stencil 
centered around jh can be listed as iVp = {(J+ l)h ■ I e J s and s = 1, 2, . . . , S} with the 
understanding that (j+l)h is a node of s th type whenever I 6 J s . This allows interaction 
between variables of multiple types. Every node (j+l )h in N^ t has a corresponding stencil 
coefficient (or weight) denoted by D ts f. Due to translational invariance, these weights 
do not change if we place the stencil at another center node j'h, hence the numbers D ts f 
do not depend on the center index j. 

We obtain the method's equation at a general node jh of the t th type by applying the 
t th stencil centered around jh to the solution values {iptj}, namely 

(38) EE%t, i+ f = 0. 

s=l le,J s 



Note that we have set all sources to zero to get a zero right hand side in (38). 

Plane waves, ip(x) = e 4fc ' 5 , are exact solutions of the Helmholtz equation with zero 
sources (and are often used to represent other solutions). Here the wave vector k is of the 
form k = u(cos(6),sin(8)) for some < 9 < 2tt representing the direction of propagation. 
The objective of dispersion analysis is to find similar solutions of the discrete homogeneous 
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system. Accordingly, we set in (38|), the ansatz 
(39) 



Vtj = Ote^ 



where k\ - Uh{cos{9), sin(#)) and a t is an arbitrary complex number associated to the t th 



variable type. We want to find such discrete wavenumbers Uh satisfying (38) 



To this end, we must solve p8j) after substituting (39) therein, namely 

(40) E^EV^ f) ^°- 



S = l UJs 



for all t = 1, 2, . . . S. Multiplying by e ik h-3h } we remove any dependence on j. Defining the 
S x S matrix F = {F ts {uih)) by 



FuM =Z D 



t,8,r 



^^(tjfjCcosfi^inS)^ )h 



UJ S 



we observe that solving (40) is equivalent to solving 

(41) detF{co h )=0. 

This is the nonlinear equation we solve to obtain the discrete wavenumber Uh correspond- 
ing to any given 9 and oj. 

5.2. Application to the DPG method. Next, we apply the above-described framework 
to the lowest order DPG stencil discussed in Section HI Since there are three different 
types of stencils (see Figure [3]), we have 5 = 3. The first type of unknowns, denoted by 



represent the DPG method's approximation to the value of <p at the nodes jh for all 



j e I 2 . The stencil of the first type is the one shown in Figure 3a The unknowns of the 
second type represent the method's approximation to the vertical components of u on 
the midpoints of horizontal edges, i.e., at all points in {KL + h/2) x KL. These unknowns 
were previously denoted by "1"' and has the stencil portrayed in Figure 3b Similarly, the 



third type of stencil and unknown are as in Figure 3c To summarize, (39) in the lowest 
order DPG case, becomes 



ip2j = Uh(xj) = a 2 e fkh '^ V% e {KL + h/2) x hZ, 

^3j = u h (x 3 ) = a 3 e fkh ^ V% e hi x {hi + h/2). 

The condensed 8x8 DPG matrices, discussed in Section |4| can be used to compute the 
stencil weights D t s j in each of the three cases, which in turn lead to the 3x3 nonlinear 



system (41) for any given propagation angle 9. 



We numerically solved the nonlinear system for u^, for various choices of 9 (propagation 
angle), r (enrichment degree), e (scaling factor in the V^-norm), and h (mesh size). The 
first important observation from our computations is that the computed wavenumbers Uh 
are complex numbers. They lie close to u in the complex plane. The small but nonzero 
imaginary parts of Uh indicate that the DPG method has dissipation errors, in addition 
to dispersion errors. The results are described in more detail below. 
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Figure 4. The curves traced out by the discrete wavevectors kh as 9 goes 
from to 7r/2. These plots were obtained using u = 1 and h = 2ir/4. 



5.3. Dependence on 9. To understand how dispersion errors vary with propagation 
angle 9, we fix the exact wave number u appearing in the Helmholtz equation to 1 (so 
the wavelength is 2ir) and examine the computed Ke(uh) for each 9. 

One way to visualize the results is through a plot of the corresponding discrete wavevec- 
tors Re{k h ) vs. k for every propagation direction 9. Due to symmetry, we only need to 
examine this plot in the region < 9 < tt/2. We present these plots for the case r = 3 
in Figure El We fix h = 27r/4. (This corresponds to four elements per wavelength if the 
propagation direction is aligned with a coordinate axis.) In Figure 4a, we plot the curve 
traced out by the endpoints of the discrete wavevectors k\- We see that as e decreases, the 
curve gets closer to the (solid) circle traced out by the exact wavevector k. This indicates 



better control of dispersive errors with decreasing e (cf. Theorem 3.1) 
In Figure 



4b 



we compare the kh obtained using the lowest order DPG method with 
the discrete wavenumbers of the standard lowest order (bilinear) finite element method 
(FEM). Clearly the wavespeeds obtained from the DPG method are closer to the exact 
u> = 1 than those obtained by bilinear FEM. However, since the lowest order DPG method 
has a larger stencil than bilinear FEM, one may argue that a better comparison is with 
methods having the same stencil size. We therefore compare the DPG method with 
two other methods which have exactly the same number of points in their stencil: (i) The 
biquadratic FEM, which after condensation has three stencils of the same size as the lowest 
order DPG method, and (ii) the conforming first order L 2 (i?) least-squares method using 
the lowest order Raviart- Thomas and Lagrange spaces (which has no interior nodes to 
condense out). While the wavespeeds from the DPG method did not compare favorably 
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(b) Dissipative errors: Plots of rj vs. e 

Figure 5. The discrepancies between exact and discrete wavenumbers as 
a function of e, when u = 1 and h = 27r/8. 



with the biquadratic FEM of (i), we found that the DPG method performs better than 
the least-squares method in (ii). 

5.4. Dependence on e and r. We have seen in Figure [4] that the discrete wavespeed 
Uh is a function of the propagation angle 9. We now examine the maximum discrepancy 
between real and imaginary parts of Uh and u over all angles. Accordingly, define 



p = max|Re(a;/ l (6 1 )) - uj\ 



r] = max |Im(u;/i ($) ) I • 



The former indicates dispersive errors while the latter indicates dissipative errors. Fixing 
u) = 1 and h = 2tt/8 (so that there are about eight elements per wavelength), we examine 
these quantities as a function of r and e in Figure |5j The first of the plots in Figures 5a 
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(a) Plot of \uJhh - wh\ for three methods 



(b) Case of DPG with r - 3 and various e 



Figure 6. Rates of convergence of \oJhh - uh\ to zero for small uh, in the 
case of propagation angle 8 = 0. 

and [5b] show that the errors decrease as e decreases from 1 to about 0.1. In view of 



Theorem |3.1[ we expected this decrease. 

However, the behavior of the method for smaller e is curious. In the remaining plots of 
Figure [5] we see that when r is odd, the errors continue to decrease for smaller e, while for 
even r, the errors start to increase as e -*■ 0. This suggests the presence of discrete effects 
due to the inexact computation of test functions. We do not yet understand it enough to 
give a theoretical explanation. 

5.5. Dependence on u. Now we examine how Uh depends on u. First, let us note that 



the matrix F in (41) only depends on uh. (This can be seen, for instance, from (37) and 
noting how the stencil weights depend on the entries of B.) Hence, we will study how 
ujhh depends on the normalized wavenumber uih, restricting ourselves to the case of 9 = 0. 



In Figure 6a, we plot (in logarithmic scale) the absolute value of Uhh-uh vs. uh for the 
standard bilinear FEM, the lowest order L 2 least-squares method (marked LS), and the 
DPG method with e = 10 _6 ,r = 3. We observe that while \uhh-uh\ appears to decrease at 
O(uh) 2 for the least squares method, it appears to decrease at the higher rate of O(uh) 3 
for the FEM and DPG cases considered in the same graph. For easy reference, we have 
also plotted lines indicating slopes corresponding to 0(uh) 2 and O(uh) 3 decrease, marked 
"quadratic" and "cubic", resp., in the same graph. 

Note that a convergence rate of \uhh-uh\ = 0(uh) 3 implies that the difference between 
discrete and exact wave speeds goes to zero at the rate 



\ujh - uj\ = us O(cohy 
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Figure 7. A comparison of discrete wavenumbers obtained by three lowest 
order methods in the case of propagation angle 9 = 0. 



This shows the presence of the so-called [2j pollution errors: For instance, as ui increases, 
even if we use finer meshes so as to maintain uh fixed, the error in wave speeds will 
continue to grow at the rate of 0(u). Our results show that pollution errors are present 



in all the three methods considered in Figure 6a The difference in convergence rates, e.g., 
whether \uh -oj\ converges to zero at the rate uj 0(uh) 2 or at the rate u 0(uh), becomes 
important, for example, when trying to answer the following question: What h should we 
use to obtain a fixed error bound for \uh - oj\ for all frequencies ul While methods with 



convergence rate uO(uh) would require h 
would only require h « ur 3 / 2 . 



to 



-2 



, methods with convergence rate uO(uhy 



Next, consider Figure pb] where we observe interesting differences in convergence rates 
within the DPG family. While the DPG method for e - 1 exhibits the same quadratic 
rate of convergence as the least-squares method, we observe that a transition to higher 
rates of convergence progressively occur as e is decreased by each order of magnitude. The 
e = 1(T 6 case shows a rate virtually indistinguishable from the cubic rate in the considered 
range. The convergence behavior of the DPG method thus seems to vary "in between" 
those of the least-squares method and the standard FEM as e is decreased. The values 
of uih considered in these plots are 2tt/2 1 for I = 1, 2, . . . , 7, which cover the numbers of 
elements per wavelength in usual practice. 

Next, we consider a wider range of uh following [21], where such a study was done for 
standard finite elements, separating the real and imaginary parts of Uhh. Our results for 
the case of 9 = are collected in Figure [7j To discuss these results, let us first recall the 
behavior of the standard bilinear finite element method (whose discrete wavenumbers are 
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also plotted in dash-dotted curve in Figure [7]). From its well-known dispersion relation 
(see e.g, [I]), we observe that if Uyji solves the dispersion relation, then 2tt -u^h also solves 



it. Accordingly, the plot in Figure 7a is symmetric about the horizontal line at height 
it. Furthermore, as already shown in [21], Uhh is real- valued in the range < uh < \/l2. 
The threshold value uh = \/l2 was called the "cut-off" frequency. (Note that in the 



regime uh > n, we have less than two elements per wavelength. Note also that v 12 > n. 



As can be seen from Figures |7a| and |7b[ in the range V 12 < uh < 6, the bilinear finite 
elements yield Uhh with a constant real part of it and nonzero imaginary parts of increasing 
magnitude. 

We observed a somewhat similar behavior for the DPG method - see the solid curves of 
Figure [7j which were obtained after calculating F explicitly using the computer algebra 
package Maple, for the lowest order DPG method, setting r = 3 and e = 0. The major 
difference between the DPG and FEM results is that u^h from the DPG method was 
not real-valued even in the regime where FEM wavenumbers were real. It seems difficult 
to define any useful analogue of the cut-off frequency in this situation. Nonetheless, we 



observe from Figures 7a and 7b that there is a segment of constant real part of value 
ir, before which the imaginary part of Uhh is relatively small. As the number of mesh 
elements per wavelength increases (i.e., as uh becomes smaller), the imaginary part of 
Uhh becomes small. We therefore expect the diffusive errors in the DPG method to be 
small when uh is small. Finally, we also conclude from Figure [7] that both the dispersive 
and dissipative errors are better behaved for the DPG method when compared to the 1? 
least-squares method. 

6. Conclusions 

We presented and analyzed the e-DPG method for the Helmholtz equation. The case 
e = l was analyzed previously in [12J. The numerical results in [12] showed that in a 
comparison of the ratio of L 2 norms of the discretization error to the best approxima- 
tion error is compared, the DPG method had superior properties. The pollution errors 
reported in [12] for a higher order DPG method were so small that its growth could not 
be determined conclusively there. In this paper, by performing a dispersion analysis on 
the DPG method for the lowest possible order, we found that the method has pollution 
errors that asymptotically grow with u at the same rate as other comparable methods. 

In addition, we found both dispersive and dissipative type of errors in the lowest order 
DPG method. The dissipative errors manifest in computed solutions as artificial damping 
of wave amplitudes (e.g., as illustrated in Figure [I]). 

Our results show that the DPG solutions have higher accuracy than an L 2 -based least- 
squares method with a stencil of identical size. However, the errors in the (lowest order) 
DPG method did not compare favorably with a standard (higher order) finite element 
method having a stencil of the same size. Whether this disadvantage can be offset by the 
other advantages of the DPG methods (such as the regularizing effect of e, and the fact 
that it yields Hermitian positive definite linear systems and good gradient approximations) 
remains to be investigated. 



DISPERSION ANALYSIS FOR DPG METHODS 21 

We provided the first theoretical justification for considering the e-modified DPG 



method. If the test space were exactly computed, then Theorem 3.1 shows that the 
errors in numerical fluxes and traces will improve as e -> 0. However, if the test space is 
inexactly computed using the enrichment degree r, then the numerical results from the 
dispersion analysis showed that errors continually decreased as e was decreased only for 
odd r. A full theoretical explanation of such discrete effects and the limiting behavior 
when e is deserves further study. 
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